;;************************************************************************
;; histfrq2.lsp 
;; contains Jan deLeeuw & Jason Bond's drawing algorithms
;;************************************************************************


(defmeth histofreq-proto :redraw ()
  (call-next-method)
  (let ((lh (+ (send self :text-ascent) (send self :text-descent)))
              (marg (send self :margin))
        (cr (send self :content-rect))
        (plot-type "Histogram")
        (right-edge (send self :canvas-width))
        (bottom-margin (- (send self :canvas-height) 
                          (fourth (send self :margin))))
        (st1)(st2)(w-st2)
        (bar-chart-legend2 (send self :bar-chart-legend2))
        )
    (when (and (not (send self :margin-flag)) (send self :n-label-lines))
          (send self :margin-flag t)
          (if (not (send self :original-bottom-margin)) 
              (send self :original-bottom-margin (fourth (send self :margin))))
          (setf (select marg 3) 
                (+ (* (- (send self :n-label-lines) 1) lh)))
          (apply #'send self :margin marg)
          (send self :margin-flag nil))
    (if (send self :bar)
        (setf plot-type (send self :legend2))
        (when (send self :poly) 
              (if (send self :prob)
                  (setf plot-type "Probability Polygon")
                  (setf plot-type "Frequency Polygon")))
        )
    (when (and (send self :histo) (not (send self :full-hist)))
          (setf plot-type (strcat "Hollow " plot-type)))
    (send self :legend2 (if bar-chart-legend2 "Bar Chart" plot-type))
    (when (and (send self :bar) (first (send self :x-axis)))
          (send self :show-cat-labels))
    (unless (send self :content-only)
            (when (not (send self :bar))
                  (setf st1 
                        (if bar-chart-legend2
                            (format nil "NBars= ~d" (ceiling (/ (send self :num-bins) 2)))
                            (format nil "NBins= ~d" (send self :num-bins))))
                  (setf st2 (if bar-chart-legend2
                                " "
                                (format nil "BinWidth= ~4,2f" (send self :binwidth))))
                  (setf w-st1 (send self :text-width st1))
                  (send self :draw-string st1 (- right-edge w-st1 2) bottom-margin)
                  (send self :draw-string st2 3 bottom-margin))
            (send self :draw-text (send self :legend1)
                  (+ (first cr) (floor (/ (third cr) 2))) 30 1 0)
            (send self :draw-text (send self :legend2)
                  (+ (first cr) (floor (/ (third cr) 2))) 45 1 0)
            (when (send self :bar)
                  (send self :draw-text (send self :variable-label 0) 
                        (floor (+ (first cr) (/ (third cr) 2))) 
                        (- (send self :canvas-height) 5) 1 0)))
    ))

(defmeth histofreq-proto :clear (&optional empty)
  (send self :start-buffering)
  (call-next-method)
  (when empty
        (send self :empty t)
        (send self :y-axis t t 0)
        (send self :x-axis t nil 0)
        (send self :variable-label 0 "")
        (send self :cat-labels (repeat " " (length (send self :cat-labels))))
        (send self :redraw)
        )
  (send self :buffer-to-screen)
  )

(defmeth histofreq-proto :show-cat-labels ()
  (let* ((x 0) 
         (yo (- (send self :canvas-height) 28  (fourth (send self :margin))));40
         (ye (+ yo 10))
         (y)
         (xy)
         (cl (send self :cat-labels))
         (cn (send self :counts))
         (tw (mapcar #'(lambda (st) (send self :text-width st)) cl))
         (wcl (max tw))
         (lh (+ (send self :text-ascent) (send self :text-descent)))
         (ncats (length cl))
         (mp (select (send self :mid-points) (* 2 (iseq ncats))))
         (mpx (mapcar #'(lambda (x) (send self :real-to-canvas x 0)) mp))
         (separation (- (first (second mpx)) (first (first mpx)) ))
         (iflag (> wcl (- (* 2 separation) 6)))
         (top-edge (select (send self :real-to-canvas 
                                 0 (max (send self :range 1))) 1))
         (bot-edge (select (send self :real-to-canvas 
                                 0 (min (send self :range 1))) 1))
         (yreal)(above-below)(fiddle)
         )
    (if (< wcl (- separation 6)) (setf ye yo))
    (mapcar #'(lambda (i)
                (if iflag 
                    (setf yreal (select cn (* 2 i)))
                    (setf yreal 0))
                (setf xy (send self :real-to-canvas (select mp i) yreal))
                (setf x (first xy))
                (if iflag
                    (setf y yo)
                    (setf y (if (evenp i) yo ye)))
                (if (not (send self :empty))
                    (send self :draw-line x (- y 3) x (- yo 7)))
                (setf y (if (evenp i) yo ye))
                (when iflag 
                      (setf above-below 0)
                      (setf fiddle 5)
                      (setf y (second xy))
                      (when (< (- (- y fiddle) (select tw i)) top-edge)
                            (setf above-below 2)
                            (setf fiddle -3)
                            (when (> (+ y (select tw i)) bot-edge)
                                  (setf above-below 0)
                                  (setf fiddle 3)
                                  (setf y bot-edge))
                            ))
                (if iflag   
                     (send self :draw-text-up (select cl i) 
                           (+ (floor (/ lh 2)) -2 x) (- y fiddle) above-below 0)
                     (send self :draw-text (select cl i) x y 1 1)))
            (iseq ncats))
    (if iflag (send self :n-label-lines 0)
        (if (= ye yo) (send self :n-label-lines 1)
            (send self :n-label-lines 2)))
    ))



;; Jan de Leeuw's drawing algorithm

(defmeth histofreq-proto :horizontal-lines ()
  (let* ((cn (send self :counts))
         (cp (send self :cut-points))
         (nn (length cn))
         (color (send self :color))
         (color? (and (send (send self :button-overlay) :color-mode) 
                      (> *color-mode* 0)))
         )
    (when (not color?) (setf color 'black))
    (dotimes (i nn)
      (let ((ci (elt cp i))
            (cj (elt cn i))
            (ck (elt cp (1+ i))))
        (when (and (send self :bar) (oddp i)) (setf ck ci))
        (send self :add-lines (list ci ck) (list cj cj) :draw nil
              :color color :width (send self :line-width))))
    ))


(defmeth histofreq-proto :poly-lines ()
  (let* ((mp (send self :mid-points))
         (cn (send self :counts))
         (nn (1- (length cn)))
         (color (send self :color))
         (color?  (and (send (send self :button-overlay) :color-mode) 
                       (> *color-mode* 0))))
    (dotimes (i nn)
      (let ((ci (elt mp i))
            (cj (elt mp (1+ i)))
            (ck (elt cn i))
            (cl (elt cn (1+ i))))
        (send self :add-lines (list ci cj) (list ck cl) :draw nil
              :color color 
              :width (send self :line-width))))))

(defmeth histofreq-proto :poly-points ()
  (let* ((mp (send self :mid-points))
         (cn (send self :counts))
         (cl (send self :cat-labels))
         (pf (if (send self :prob) "P" "f"))
         (j 0)
         (color 'red)
         (color?  (and (send (send self :button-overlay) :color-mode) 
                       (> *color-mode* 0))))
    (send self :add-points mp cn :draw nil :color color :symbol 'disk)
    (mapcar #'(lambda (i)
                (when (oddp i) (setf j (1+ j)))
                (send self :point-label i
                      (if (send self :bar)
                          (format nil "~a=~d"pf (fuzz (elt cn i)))
                          (format nil "~a=~d midpt=~3,2f" pf 
                                  (fuzz (elt cn i)) (elt mp i))))
                (when (and (send self :bar) (oddp i)) 
                      (send self :point-state i 'invisible)))
            (iseq (length mp)))
            ))

(defmeth histofreq-proto :hollow-vertical-lines ()
  (let* ((cn (send self :counts))
         (cp (send self :cut-points))
         (nn (length cp))
         (color  (send self :color))
         (color? (and (send (send self :button-overlay) :color-mode) 
                      (> *color-mode* 0))))
    (send self :add-lines (repeat (first cp) 2)
          (list 0 (first cn)) :draw nil
          :color color 
          :width (send self :line-width))
    (send self :add-lines (repeat (first (last cp)) 2)
          (list 0 (first (last cn))) :draw nil
          :color color 
          :width (send self :line-width))
    (dolist (i (1+ (iseq (- nn 2))))
      (let ((ci (elt cp i))
            (cj (elt cn (1- i)))
            (ck (elt cn i)))
        (send self :add-lines (list ci ci)
              (list (min cj ck) (max cj ck)) :draw nil
              :color color 
              :width (send self :line-width))))))

(defmeth histofreq-proto :full-vertical-lines ()
  (let* ((cn (send self :counts))
         (cp (send self :cut-points))
         (nn (length cp))
         (color  (send self :color))
         (color? (and (send (send self :button-overlay) :color-mode) 
                       (> *color-mode* 0))))
    (dolist (i (1+ (iseq (- nn 2))))
      (let ((ci (elt cp i))
            (cj (elt cn (1- i)))
            (ck (elt cn i)))
        (send self :add-lines (list ci ci)
              (list 0 (min cj ck)) :draw nil
              :color color 
              :width (send self :line-width))))))



;;  Jan de Leeuw's Compute Methods

(defmeth histofreq-proto :recompute ()
  (let ((counts (send self :make-counts)))
    (when (> (length (send self :data)) (length counts) 2)
          (send self :counts (send self :make-counts))
          (send self :cut-points (send self :make-cut-points))
          (send self :mid-points (send self :make-mid-points)))
    (> (length (send self :data)) (length counts) 2)))

(defmeth histofreq-proto :make-counts ()
  (let ((data (send self :data))
        (binwidth (send self :binwidth))
        (origin (send self :origin)))
    (send self :count-in-bins data binwidth origin)))

(defmeth histofreq-proto :make-cut-points ()
  (let* ((data (send self :data))
         (binwidth (send self :binwidth))
         (origin (send self :origin))
         (cn (length (send self :counts)))
         (minx (floor (/ (- (min data) origin) binwidth)))
         (maxx (ceiling (/ (- (max data) origin) binwidth)))
         (seq (if (= (1+ cn) (- (1+ maxx) minx))
                  (iseq minx maxx)
                  (if (= cn (- (1+ maxx) minx))
                      (iseq minx (1+ maxx))
                      (iseq minx (1- maxx))))))
    (+ origin (* binwidth seq))))

(defmeth histofreq-proto :make-mid-points ()
  (let ((cp (send self :cut-points)))
    (/ (+ (butlast cp) (rest cp)) 2)))

;;  Initial Constructors

(defmeth histofreq-proto :make-origin ()
  (send self :origin 
        (mean (send self :data))))

(defmeth histofreq-proto :make-binwidth ()
  (send self :binwidth
        (cond 
          ((send self :histo)
           (send self :opt-histo-binwidth))
          ((send self :poly)
           (send self :opt-poly-binwidth))
          ((send self :bar) 1)
          (t (error "Unknown plot type"))))
  (send self :binwidth (/ (round (* 10 (send self :binwidth))) 10))
  )

(defmeth histofreq-proto :opt-histo-binwidth ()
"See Scott, Biometrika, 66, 1979, 605-610"
  (let (
        (data (send self :data))
        )
    (* (standard-deviation data) 
       (^ (* 24 (sqrt pi) (/ (length data))) (/ 3)))
    ))

(defmeth histofreq-proto :opt-poly-binwidth ()
"See Scott, JASA, 80, 1985, 348-454"
  (let (
        (data (send self :data))
        )
    (* 2 (standard-deviation data) 
       (^ (* (/ 40 49) (sqrt pi) 
             (/ (length data))) (/ 5)))))
    

(defmeth histofreq-proto :count-in-bins (data binwidth origin)
  (if (listp binwidth)
      (progn 
       (let (
             (beg-bin (butlast binwidth))
             (end-bin (rest binwidth))
             (bcnt (repeat 0 (1- (length binwidth))))
             )
         (mapcar #'(lambda (x) 
                     (let ((posn (first (which (if-else (< beg-bin x end-bin) 
                                                        t nil)))))
                       (setf (elt bcnt posn) (1+ (elt bcnt posn))))) data)
         bcnt))

      (progn (let* (
                    (indx (floor (/ data binwidth)))
                    (jndx (- indx (min indx)))
                    (nbin (1+ (max jndx)))
                    (bcnt (repeat 0 nbin))
                    )
               (mapcar #'(lambda (x) (setf (elt bcnt x) (1+ (elt bcnt x)))) jndx)
               bcnt
               ))))

(defmeth histofreq-proto :spacing (data)
  (min (abs (difference (sort-data data)))))

(defmeth histofreq-proto :num-bins ()
  (length (send self :counts))
  )

(defmeth histofreq-proto :num-cuts ()
  (length (send self :cut-points))
  )


;; Forrest Young's Normal and Kernel Density curve fitting.
;; Not least squares. Simply adds curves with the data's mean and stdv.

(defmeth histofreq-proto :choose-density ()
  (when (> (send self :num-points) 2)
        (cond 
          ((not (send self :slot-value 'dens-dialog))
           (let* ((den-fun-state (list nil nil nil))
                  (monopoly-arguments nil)
                  (title  (send text-item-proto :new "Choose Curve"))
                  (Kernel-text (send text-item-proto :new "Kernel:"))
                  (normal (send toggle-item-proto 
                    :new "Normal Density"
                    :value (select den-fun-state 0)
                    :action #'(lambda () 
                                (send self :switch-add-normal))))
                  (kernel-type (send choice-item-proto
                    :new(list "Bisquare" "Gaussian" "Triangular" "Uniform")
                    :value 0
                    :action #'(lambda () 
                                (send self :put-kernel-type))))
                  (kernel (send toggle-item-proto 
                    :new "Kernel Density"
                    :value (select den-fun-state 1)
                    :action #'(lambda ()
                      (send self :kernel-type (send kernel-type :value))
                      (send self :switch-add-kernel))))            
                  (monopoly (send toggle-item-proto 
                    :new "MonoPoly Density"
                    :value (select den-fun-state 2)))
                  (dialog (send dialog-proto 
                    :new (list title normal 
                               (list (list kernel 
                                           (list kernel-text
                                                 kernel-type))))
                    :title "Distribution Curve"))
                 )
             (defmeth self :put-kernel-type ()
               (send self :kernel-type (send kernel-type :value))
               (when (send self :show-kernel)
                     (send self :switch-add-kernel)
                     (send self :switch-add-kernel)))
             (send self :add-subordinate dialog)
             (send self :slot-value 'dens-dialog dialog)))
          (t
           (send (send self :slot-value 'dens-dialog) :show-window)))))
          

(defmeth histofreq-proto :add-normal (&key (color 'red) line-width (draw t))
        (let* (
               (current-variable (first (send self :content-variables)))
               (ndim (- (send self :num-variables) 1))
               (var (send self :data))
               (mu (mean var))
               (s (standard-deviation var))
               (range (send self :range current-variable))
               (x (rseq (first range) (second range) 50))
               (y (normal-dens (/ (- x mu) s)))
               (freq (send self :point-coordinate  1
                           (iseq (send self :num-points))))
               (y (* y (/ (max freq) (max y))))
               (overlay (first (send self :slot-value 'overlays)))
               )
          (when (not line-width) (setf line-width *line-width*))
          (send self :add-lines (append (repeat (list x) ndim) (list y))
                :draw nil
                :color color ;(send self :normal-curve-color)
                :width line-width)
          (when draw (send self :redraw))))


(defmeth histofreq-proto :add-kernel 
                         (type-number &key (color 'green) line-width (draw t))
  (let* ((type (case type-number
                 (0 'b)
                 (1 'g)
                 (2 't)
                 (3 'u)))
         (current-variable (first (send self :content-variables)))
         (npts (send self :num-points))
         (ndim (- (send self :num-variables) 1))
        ;(varx (send self :point-coordinate current-variable (iseq npts)))
         (varx (send self :data))
         (x (rseq (min varx) (max varx) 50))
         (maxyaxis (max (send self :range 1)))
         (y (second (kernel-dens varx :xvals 50 :type type)))
         (maxy (max y))
         (x-at-maxy (select x (first (last (order y)))))
         (freq (send self :point-coordinate  1
                           (iseq (send self :num-points))))
         (y (* y (/ (max freq) (max y))))
         (overlay (first (send self :slot-value 'overlays)))
         )
    (when (not line-width) (setf line-width *line-width*))
    (send self :add-lines (append (repeat (list x) ndim) (list y)) 
          :draw nil
          :color color ;(send self :kernel-curve-color)
          :width line-width)
    (when draw (send self :redraw))))

(defmeth histofreq-proto :switch-add-normal 
      (&key (color 'red) (kcolor 'green) line-width (draw t))
      (send self :show-normal (not (send self :show-normal)))
      (cond
        ((send self :show-normal)
         (send self :add-normal
               :draw draw 
               :color 'red ;(send self :normal-curve-color)
               :line-width line-width))
        (t
         (send self :start-buffering)
         (send self :clear-lines)
         (send self :make-plot)
         (when (send self :show-kernel) 
               (send self :add-kernel 
                    (send self :kernel-type) 
                     :draw draw
                     :color 'green ;(send self :kernel-curve-color)
                     :line-width line-width))
         (send self :buffer-to-screen))))
    
(defmeth histofreq-proto :switch-add-kernel (&key (color 'red) line-width (draw t))
      (send self :show-kernel (not (send self :show-kernel)))
      (cond 
        ((send self :show-kernel)
         (send self :add-kernel 
              (send self :kernel-type) 
               :draw draw
               :color 'green ;(send self :kernel-curve-color) 
               :line-width line-width))
        (t
         (send self :start-buffering)
         (send self :clear-lines)
         (send self :make-plot)
         (when (send self :show-normal) 
               (send self :add-normal 
                     :draw draw
                     :color 'red ;(send self :normal-curve-color) 
                     :line-width line-width))
         (send self :buffer-to-screen)
         )))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Fit Normal Curve to HistoFreq -- written by Jan de Leeuw
;; Divide-by-zero problems fixed by Forrest Young
;;
;; Inspired by Brown and Hwang, American Statistician, 1993,
;; 251-255, see also, 1994, 353-354. It is far more general, however.
;;
;; Takes a scatterplot-proto, extracts all the linestarts, then
;; reconstructs from this all the line segments in the plot, eliminating
;; the vertical segments in the process. These segments can result from
;; a (normalized) histogram, or from a function plot, or just from
;; an arbitrary bunch of linestarts. Some lines can be horizontal,
;; some can run from right to left. Segments are coded by their
;; endpoints on the horizontal axes, intercepts, and slopes.
;;
;; The program then minimizes
;;
;; sum_{i=1}^n int_{x_i}^{y_i}(a_i+b_i z-cf(z,m,s))^2dz
;;
;; where x_i,y_i are the endpoints of segment i, a_i is the
;; intercept of the correspoding line and b_i is its slope,
;; and f(z,m,s) is the density of a normal with mean m and
;; standard-deviation s.
;;
;; We minimize either over m,s,c (default), i.e. we also
;; scale the fitted density, or we minimize over m,s fixing
;; c at 1. This can be used to fit a normal to a histogram,
;; but also to fit a normal to the plot of a gamma density,
;; and so on.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
          

(defmeth histofreq-proto :refit-normal ()
  (let* ((aa (fit-normal (send self :convert-linestarts)))
         (nn (send self :num-lines))
         (ll (+ (- nn 50) (iseq 50))))
    (send self :linestart-coordinate 0 ll (first aa))
    (send self :linestart-coordinate 1 ll (second aa))
    (if (send self :normal)
        (send self :toggle-normal-plot))))

(defmeth histofreq-proto :toggle-normal-plot ()
  (let* ((nn (send self :num-lines))
         (ll (+ (- nn 50) (iseq 49))))
    (if (send self :normal)
        (send self :linestart-next ll (1+ ll))
        (send self :linestart-next ll nil))
    (send self :adjust-to-data)))

(defmeth histofreq-proto :normal-lines ()
  (let ((ll (fit-normal (send self :convert-linestarts))))
    (when ll
          (send self :add-lines (first ll) (second ll)
                :draw nil))
    ll))
#|
(defun fit-normal
  (lst &key (mean 0) (sigma 1) (scale t) (draw t) (verbose nil))
"Args: x
Given a list of linestarts, which define a number of
segments, this fits a normal curve to the segments by minimizing mean
square error."
  (let* ((arg (list mean sigma))
         )
(break)
    (setf arg 
          (newtonmax #'(lambda (arg) (nfl-func lst arg scale))
                     arg :verbose verbose))
    (if verbose  (format t "~,10f ~%" (nfl-func lst arg scale)))
    (let* ((ccc (normal-fit-collect-coefs lst (first arg) (second arg)))
           (c11 (second ccc))
           (c22 (third ccc))
           (srt (mapcar #'first lst))
           (end (mapcar #'second lst))
           (xes (rseq (min srt) (max end) 50)))
      (if (= c22 0) 
          (format t "~%Optimization method cannot fit Normal Curve")
          (list xes (mapcar
                     #'(lambda (x)
                         (let ((fac (if scale (/ c11 c22) 1))
                               (mn (first arg))
                               (sg (second arg)))
                           (* fac (/ (normal-dens (/ (- x mn) sg)) sg))
                           )) xes))))))

(defun nfl-func (lst arg-list scale)
  (let ((result (normal-fit-loss-function
                 lst (first arg-list) (second arg-list) scale))
        )
    (when result (setf result (- result)))))

(defmeth scatterplot-proto :convert-linestarts ()
"Args: lines
The list of linestarts is converted to a list of endpoints,
slopes, and intercepts."
(let* ((nn (send self :num-lines))
       (sg nil))
  (dotimes (i nn)
    (let ((j (send self :linestart-next i)))
      (if j (let ((x0 (send self :linestart-coordinate 0 i))
                  (y0 (send self :linestart-coordinate 1 i))
                  (x1 (send self :linestart-coordinate 0 j))
                  (y1 (send self :linestart-coordinate 1 j)))
              (if (/= x0 x1)
                  (let ((a (/ (- (* y0 x1) (* y1 x0))
                              (- x1 x0)))
                        (b (/ (- y1 y0) (- x1 x0))))
                    (setf sg (append sg (list x0 x1 a b)))))))))
  (split-list sg 4)
  ))

(defun normal-fit-collect-coefs (lst mn sg)
  (let* ((srt (mapcar #'first lst))
         (end (mapcar #'second lst))
         (int (mapcar #'third lst))
         (slp (mapcar #'fourth lst))
         (c00 (sum (+ (* (^ slp 2) (/ (- (^ end 3) (^ srt 3)) 3))
                      (* slp int (- (^ end 2) (^ srt 2)))
                      (* (^ int 2) (- end srt)))))
         (df1 (/ (- srt mn) sg))
         (df2 (/ (- end mn) sg))
         (c11 (sum (- (* (+ int (* mn slp))
                         (- (normal-cdf df2)
                            (normal-cdf df1)))
                      (* sg slp
                         (- (normal-dens df2)
                            (normal-dens df1))))))
         (c22 (sum (/ (- (normal-cdf (* df2 (sqrt 2)))
                         (normal-cdf (* df1 (sqrt 2))))
                      (* 2 sg (sqrt pi))))))
    (list c00 c11 c22)
    ))

(defun normal-fit-loss-function (lst mn sg sc)
  (let* ((ccc (normal-fit-collect-coefs lst mn sg))
         (c00 (first ccc))
         (c11 (second ccc))
         (c22 (third ccc)))
    (cond
      ((= 0 c22) 
    ;  (format t "~%Normal-fit-loss-function: Avoiding division by zero")
       0)
      (t
       (if sc (- c00 (/ (^ c11 2) c22))
           (- (+ c00 c22) (+ c11 c11)))))))


|#